Mode space approach for tight-binding transport simulations in graphene nanoribbon 
field-effect transistors including phonon scattering 
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In this paper, we present a mode space method for atomistic non-equilibrium Green's function 
simulations of armchair graphene nanoribbon FETs that includes electron-phonon scattering. With 
reference to both conventional and tunnel FET structures, we show that, in the ideal case of a 
smooth electrostatic potential, the modes can be decoupled in different groups without any loss of 
accuracy. Thus, inter-subband scattering due to electron-phonon interactions is properly accounted 
for, while the overall simulation time considerably improves with respect to real-space, with a 
speed-up factor of 40 for a 1.5-nm-wide device. Such factor increases with the square of the device 
width. We also discuss the accuracy of two commonly used approximations of the scattering self- 
energies: the neglect of the off-diagonal entries in the mode-space expressions and the neglect of the 
Hermitian part of the retarded self-energy. While the latter is an acceptable approximation in most 
bias conditions, the former is somewhat inaccurate when the device is in the off-state and optical 
phonon scattering is essential in determining the current via band-to-band tunneling. Finally, we 
show that, in the presence of a disordered potential, a coupled mode space approach is necessary, 
but the results are still accurate compared to the real-space solution. 



I. INTRODUCTION 

Graphene nanoribbons (GNRs) have been proposed in 
recent years as a possible replacement for silicon in the 
future generation of field effect transistors Despite 
the outstanding challenges in producing GNRs with con- 
trolled width and edges and the resulting degradation of 
the graphene intrinsic mobility, devices with large on-off 
current ratio have been successfully demonstrated Q. 

Theoretical performance of GNRFETs have been 
widely investigated by numerous simulation studies (e.g 
Refs. The state-of-the-art approach for modeling 

the electronic properties of GNRs is based on an atom- 
istic tight-binding (TB) Hamiltonian with a pz orbital 
basis set. The transport problem is usually solved within 
the nonequilibrium Green's function formalism (NEGF) 
[lo| . which provides a rigorous framework for including 
incoherent scattering processes in the quantum descrip- 
tion. TB simulations of GNRFETs including phonon 
scattering have been reported too [TT| - [l3 | . Those simula- 
tions were performed using a real space (RS) approach. 
On the other hand, more efficient mode space (MS) meth- 
ods would be preferable for use in intensive device simu- 
lations. 

The MS approach is well established with reference to 
an effective mass (EM) Hamiltonian [l^ [l^ . It is based 
on the expansion of the Green's functions in terms of 
the transverse eigenfunctions (modes) and it is most ef- 
ficient when quantum confinement is relatively strong in 
the transverse plane, so that only few of the lowest sub- 
bands are occupied. The so-called coupled mode space 
(CMS) is the general method, while a more efficient ver- 
sion, the uncoupled mode space (UMS), can be adopted 
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when the transverse potential profile is slightly varying 
along the transport direction. The inclusion of electron- 
phonon scattering in the MS EM model is also well known 
[17,-^1. 



On the other hand, the MS approach is not gener- 
ally applicable to a TB Hamiltonian. Due to the non- 
separable Hamiltonian, even when the potential is uni- 
form along the longitudinal direction (so that the longi- 
tudinal wavevector k is conserved), the wavefunctions at 
different k in the same subband are in general different. 
In other words, by defining the modes as the transverse 
part of the wavefunctions at a particular fc, the wave- 
functions at the generic k are a linear combination of 
them. This is similar to what occurs in the case of a fc • p 
Hamiltonian [g^]. Nevertheless, for the particular case of 
the TB Hamiltonian of GNRs with armchair edges, we 
have previously shown by numerical calculation that a 
MS method is actuall y p ossible, since each subband con- 
tains only few modes [23j . A formal derivation was given 
in Ref. 2^ using analytically defined modes. 



In this paper, we propose a novel MS method for the 
inclusion of graphene acoustic phonon (AP) and optical 
phonon (OP) scattering within the NEGF formalism, and 
test its accuracy with respect to RS in both the cases of 
smooth and disordered potentials. Two approximations 
commonly found in the literature for treating the scatter- 
ing self-energy are also evaluated. The paper is organized 
as follows. Sec. HI] reviews the standard RS TB method 
as well as the MS TB method from Ref. [2^ modified so 
as to include phonon scattering. Simulation results of 
GNRFETs are presented in Sec. IIIII and conclusions are 
finally drawn in Sec. IIVI 
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FIG. 1. One-dimensional elementary cell of an Na = 13 arm- 
chair GNR. 



II. MATHEMATICAL MODEL 
A. Real space formulation 

The atomic structure of an armchair GNR unit cell is 
shown for reference in Fig. [TJ The GNR width is denoted 
by the number of dimer lines Na- We adopt the TB 
Hamiltonian model proposed in Ref. [251, which is based 
on a set of orthogonal Pz orbitals, one for each atom, 
with hopping integrals limited to first nearest-neighbor 
orbitals and of value t for internal atom pairs, and t{l + 6) 
for atom pairs located along the edges of the GNR {t — 
—2.7 eV and 6 = 0.12). The onsite energies are set equal 
to the value of the electrostatic potential energy at each 
atomic site. Unless stated otherwise, the electrostatic 
potential is calculated by self-consistently solving the 3D 
Poisson equation. 

In order to calculate electron and hole densities within 
the NEGF formalism, the following equations need to be 
solved for the retarded, the lesser and the greater Green's 
function matrices G^, and G^, respectively, for each 
energy E 

[{E + iO+)I-H'^-T,^{E)]G^{E)=I (1) 



G<{E) = G^iE)T,<{E)G^{E) 



G>{E) = G"{E)T,>{E)G'^{E) 



(2) 



(3) 



where H"^ denotes the restriction of the Hamiltonian ma- 
trix H to the GNR portion inside the simulation domain 
and G^ = G^l Besides, S^(£;) = + Sf + 

Sg(£') and similarly for 1:>{E) and where S^y], 

are the self-energies for the semi-infinite source and drain 
leads, which are assumed in thermodynamic equilibrium, 
while Sp '' account for phonon scattering. Calculation 
of both ([2]) and ([3]) is actually not needed, since either 
one of the two equations can be replaced by the identity 
G> - G< = G« - G^. 

Phonon scattering is treated within the self-consistent 
Born approximation. The model accounts for both OP 
and AP scattering, the latter within the elastic and high 
temperature limit, resulting in the following expressions 



of the lesser and greater phonon self energies 
S<(S) = [DapG<iE) + DopiNop + 1)G<{E + hiOop) + 
+ DopNopG<{E - fii^op)] o I (4) 



S>(S) = [DapG>iE) 



-^op op 



DopNopG' 



1)G>{E -huoop) 



'>{E + nujop)\oI (5) 

where o indicates element-by-element matrix multiplica- 
tion, Nop is the OP occupation number, and Dap and 



Dop are given by 



_DlkBT 

^ap — 



D. 



op 



SrrirUJn, 



(6) 



with Dac = 16 eV the AP deformation potential, Vg — 
2 X 10^ cm/s the sound velocity in graphene, T the tem- 
perature, rric the carbon atomic mass. Do = 10^ eV/cm 
and hujop — 160 meV the zone-boundary OP deformation 
potential and energy, respectively. The retarded phonon 
self-energy is computed using the identity 



4-00 



s«(£;) 



i f ^>{E')-^<{E': 



2tt 



E + iO+ - E' 



dE' 



+ 00 



- 2 2^{E-E') '^^ 

— oo 

where the symbol P stands for the principal part of the 
integral. Unless stated otherwise, we include the con- 
tribution to the principal part integral (Hermitian part 
of Sp) due to AP scattering, while we neglect the one 
due to OP scattering. See Appendix 1X1 for details on the 
numerical implementation of the principal part integral. 

The equations above, namely ([I])-© together with (jl]), 
([5]) and ([7]), are iteratively solved. It is worth noticing 
that if* and S^ - ) can be ordered in a block tridiago- 
nal and block diagonal form, respectively, each diagonal 
block having the size of the number of atoms in a spe- 
cific transverse layer, or slab, of the GNR. Hence, the 
Recursive Green Function algorithm (RGF) [2^ [23| can 
be used to solve ©-(IS]). In this paper, we consider each 
slab to be made of four rows of atoms, i.e. equal to the 
GNR unit ceh (Fig.[T]). 



B. Mode space formulation 

In this section the MS approach presented in Ref. [23l . 
and in that paper discussed with reference to coherent 
transport simulations of GNRs, is extended to include 
phonon scattering. 

The MS formulation starts by defining a set of or- 
thonormal vectors, or modes, 4>{i) for each slab i of the 
device. We choose the modes of the generic slab as the 
eigenvectors computed at fc = of a fictitious GNR ob- 
tained by the periodic repetition (potential energy in- 
cluded) of that slab along the longitudinal direction. In 
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the absence of topological differences between the slabs 
(such as edge irregularities or internal vacancies), the 
modes are the solution of the eigenvalue problem 



H 



■i+i 



(f,"''{i)^e"'(p"'{i) (8) 



where Hi^j is the Hamiltonian block between slab i and 
J. Denoting by v{i) the matrix whose columns are the 
modes of layer z, i.e. v{i) — [4>^{i) ■ ■ ■ cf)"^{i) ■ ■ ■], a block- 
diagonal transformation matrix V satisfying = I 
is constructed 



V = 



(v{l) 
v{2) 

V ■•• 



\ 



(9) 



with Ng equal to the number of slabs. 

In the following, we indicate with a tilde the MS matri- 
ces to distinguish them from the RS matrices. The CMS 
method consists in approximating (and similarly for 
G< and G>) as 



G^{E) ~ VG^iE)V^ 
where G^ is the solution of 

{E + iO+)J H'^- T,^{E) \ G'^iE) = I 



with 



(10) 

(11) 

(12) 
(13) 



Eq. (jlip is the MS version of ([T]). Analogous consider- 
ations apply to the other equations. Eq. (fTO|) would be 
exact if V was a square matrix. In practice, a mode trun- 
cation is performed so that H"^ and have a smaller 
size than their corresponding RS matrices and the solu- 
tion of ([TT]) instead of ([T]) is computationally advanta- 
geous. The computational time of the RGF algorithm 
scales as 0{NyNs), where Ny is the matrix block size 
[271] . which is equal to the slab size for RS {2Na accord- 
ing to our choice of the slab size) and to the number Nj^ 
of selected modes per slab for CMS. Thus, the speed- 
up of CMS compared to RS is a factor of the order of 

An algorithm to select the modes to retain in the V 
matrix was presented in Ref. :23. It is based on identify- 
ing, among the modes calculated with zero electrostatic 
potential, the minimum set of modes that allow to repro- 
duce with sufficient accuracy the subbands that lie in the 
energy range of interest. The mode indexes so identified 
are then used to select the actual modes calculated with 
the non-null electrostatic potential. As shown in Fig. [5] 
for a Na — 13 GNR, it turns out that each one of the 
lowest conduction subbands at zero potential can be well 
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FIG. 2. Subband structure of a iVa = 13 GNR computed 
with RS (dashed lines) and MS (solid lines) using two difler- 
ent groups of modes (left and right). The modes included in 
each group are the ones that correspond to the eigenvalues at 
fc = indicated with circles in each figure. The electrostatic 
potential is set to zero. Sa^c is the length of the GNR unit 
cell, with flee the carbon interatomic distance. 



reproduced by just four modes, which correspond to the 
eigenvalues at A: = belonging to (i) the considered sub- 
band, (ii) the valence subband symmetrical to it, and 
(iii-iv) their respective folded continuations. Symmetri- 
cal considerations apply to the highest valence subbands. 
GNRs with different Na behave similarly. Since the sets 
of four modes described above, hereafter referred to as 
groups of modes, are disjoint from each other, a more 
efficient MS method, here called uncoupled group mode 
space (UGMS), can be obtained from p0)) - p3)) by ne- 
glecting the coupling between modes belonging to dif- 
ferent groups. Let Vh^i) be the matrix whose columns 
are the modes in layer i of group b only. A transforma- 
tion matrix Vf, similar to ([9]) is constructed by using the 
matrices vi, as diagonal blocks. The UGMS method is 
derived by further approximating (flU)) as 



G^{E)^yV,G^{E)V,'^ 



6=1 



(14) 



where Ng is the number of the considered groups and G^ 
is the solution of 



\e + iO+)/ - Ht S«(i?)l G^{E) = I (15) 



with 



sf (s) - yfctsfl(^)v, 



(16) 
(17) 



In the UGMS method, the computational cost of the 
RGF algorithm scales as 0{NgN^Ns), where Ny = A 



4 



is the number of modes in each group. Hence, if the 
same number of modes A^^^ = is used, the speed-up 
of UGMS compared to CMS scales as iV^ , which can be 
sizeable for wide ribbons, since Ng is proportional to the 
number of subbands that contribute to transport, which, 
in turn, is proportional to the GNR width. 

In the following, we refer to the UGMS formulation, 
since the CMS one can be recovered from Eqs. (|14I) - PT)) 
by considering all the selected modes as belonging to the 
same group and by setting Ng = 1. Note that and 
can be directly computed in MS without making use 

of (fT7|) . Instead, the calculation of Sp '' is more com- 
plicated. In particular, in the MS representation these 
matrices are no longer diagonal, but just block diagonal. 
For example, as far as Sp is concerned, from the equa- 
tion for S< analogous to (|17p . using (j4]) and considering 
for simplicity the contribution of APs only (first term at 
the right-hand-side), with G< replaced by the equation 
analogous to (|14p . the following MS expression for the el- 
ement relative to modes m and m' of the diagonal matrix 
block associated with slab i and group b is derived 

^b,m7n'{i,i; E) = Dap ^ Fb.mm' {"i) Gb' ,nn' {i, i', E) 
b' .n,n' 

(18) 

where the symbols < and P have been dropped for 
brevity and the form factor F is defined by 

a 

(19) 

(n and n' are mode indexes within group b' , a the atom 
index within the slab). Identical considerations apply 
to the OP terms and to the other types of phonon self- 
energies. It should be noticed that all the considered 
modes give their contribution to Eb^mm', not only those 
belonging to the same group b of modes m,m' . This 
means that intersubband scattering between all consid- 
ered modes is taken into account, without the necessity 
of going to full RS simulations. It must also be noticed 
that in the solution process the different mode groups 
can be treated independently, which is a great advantage 
compared to RS and also to CMS: at every iteration pass, 
once the Green functions for all groups have been calcu- 
lated, the self-energies are updated using the equations 
like (ITS)) collecting the contributions from all modes. 

In the next section the MS and RS solutions will be 
compared in the presence of phonon scattering. A fur- 
ther approximation will be tested, which consists in sim- 
plifying the form factor as 

-^&,mm'(*) - ^yn,m'5n,n' ^ \vb,am{i)\^ W ,an{i)\^ (20) 
a 

where 5 is the Kronecker delta, so that only the diagonal 
entries are kept in (|18p and similar equations. This ap- 
proximation is usually adopted in the context of the CMS 
EM approach due to its higher efficiency [13, El, HI • The 
results obtained with expressions and (PU)) will be 
named "MS2" and "MSI", respectively. 



III. RESULTS 

All simulated devices have a double-gate structure with 
Si02 1-nm-thick top and bottom oxides. As mentioned 
before, the source and drain are assumed to be semi- 
infinite leads. The portions of the source and drain re- 
gions included in the simulation domain are 10 nm long 
and uniformly doped. The channel is assumed to be in- 
trinsic. Unless noted otherwise, the following device pa- 
rameters are used: Na = 13 (corresponding to a ribbon 
width ~ 1.5 nm), gate length Lq = 17 nm, and dop- 
ing concentration in the source and drain regions equal 
to 10^^ dopants per carbon atom. Also, unless stated 
otherwise, the MS method simulations are performed us- 
ing the uncoupled group approximation (UGMS). For the 
reference Na = 13 GNR, the 2 groups of 4 modes of Fig. [2] 
are used (8 modes out of a total of 26 that correspond to 
the full RS solution) . Both conventional FET and tunnel 
FET (TFET) devices are simulated. The latter have re- 
ceived great attention in recent years for their potential 
in low-power applications [s^]. The two structures only 
differ in the type of doping of the source region. While 
the drain is always n-type, the source is n-type and p- 
type for conventional FETs and TFETs, respectively. 

We consider only ideal GNRs with perfect edges and 
no internal vacancies. Nevertheless, in the following, we 
also study the effect of adding a disorder potential to 
the TB Hamiltonian of the ideal GNR. Such disorder 
potential can, in principle, mimic the perturbation effect 
of impurities or the substrate. 



A. Smooth electrostatic potential 

The case of no disorder potential is treated first. 

The / vs. Vgs ("turn-on") characteristics of a con- 
ventional n-i-n FET at Vds — 0.1 and 0.8 V are plotted 
in Fig. [3l They are computed with the RS and the MS 
methods. For each method, we consider three transport 
conditions: (i) without scattering (i.e. ballistic trans- 
port), {a) in the presence of only AP scattering, and 
[iii) in the presence of both AP and OP scattering. In 
the case of phonon scattering and MS approach, the sim- 
plified expression of the form- factor in (|20p is used (MSI). 

As regards the effect of phonon scattering and with ref- 
erence to the RS results, it can be seen that AP scattering 
has only a limited effect on the current at this channel 
length, resulting in a ballisticity ratio (i.e. ratio between 
current in the presence of phonon scattering and ballis- 
tic current) of 0.8 at Vgs = 0.8 V for both Vds values. 
When also OP scattering is included, the current at high 
Vgs, i-G. on-state current, is only slightly decreased for 
the largest Vds value. Similar findings were reported in 
Ref. [12. On the other hand, OP scattering is responsi- 
ble for an increase of the minimum off-state current by 
a few orders of magnitudes when Vds is low. This effect 
is caused by energy relaxation through emission or ab- 
sorption of optical phonons, which favors band-to-band 
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N =13, L„=17nm 

a ' G 




FIG. 3. Turn-on characteristics of a n-i-n FET with W — 
1.5 nm at Vds = 0.1 V (top) and Vds = 0.8 V (bottom). Re- 
sults are obtained by the RS and UGMSl methods for differ- 
ent transport conditions: ballistic transport, AP scattering, 
and both AP and OP scattering. 
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FIG. 4. Current spectrum Jb{x,E) obtained by UGMSl 
method for mode group b = 1 (left) and h = 2 (right) for the 
device in Fig. [3] with AP scattering and at the bias indicated 
in the figure. Blue (red) color means low (high) density. The 
white solid (dashed) line is the profile of the first (second) con- 
duction subband. The source Fermi level is at i5 = 0. The 
inset shows the integrals over energy Ih(x) = / Jb{x, E)dE 
and the conservation of the total current Ii{x) + l2{x). 



tunneling (BTBT) and shifts toward positive Vqs values 
the onset of the ambipolar conduction by BTBT, similar 
to what occurs in carbon nanotubes j28| . 

With regard to the MS results, we do not observe any 
significant discrepancy with respect to RS, except for the 
curve with OP scattering at Vds = 0.1 V, close to the 
Vcs point of minimum current. This discrepancy will be 
discussed later in the text. 

In order to test the validity of the MS method in de- 




X (nm) X (nm) 



FIG. 5. Same as in Fig.[4]but with both AP and OP scattering 
and at Vas = -0.2 V and Vds = 0.1 V. The white solid 
(dashed) lines are the profiles of the first (second) pairs of 
conduction and valence subbands. 



scribing inter-subband scattering, we separately plot in 
Fig. |4] the MS current spectrum for the first and second 
group of modes (corresponding to the first and second 
subband, respectively) at Vgs — 0.8 V and Vds = 0.1 V, 
in the case with only AP scattering. It can be seen that 
both subbands contribute to current. In addition, despite 
the absence of inelastic scattering processes in the simu- 
lation, the distribution of current over energy of the sin- 
gle subband is not conserved when moving from source to 
drain, indicating that some of the carriers are transferred 
from one subband to the other. Indeed, thanks to phonon 
scattering, part of the electrons injected from the source 
in the second subband reach the drain by traversing the 
channel in the first subband, where the energy barrier is 
lower. The total current is properly conserved as shown 
in the inset of Fig. 2] Inter-subband scattering and the 
conservation of the total current are also evident from 
simulations including OP scattering: see Fig. [5l which 
corresponds to a bias point where transport is dominated 
by phonon-assisted BTBT. 

Next, we consider the TFET architecture. Fig. [6]-top 
shows the / vs. Vgs characteristics at Vds = 0.4 V 
for the p-i-n counterpart of the device in Fig. [31 calcu- 
lated with the different methods and by including differ- 
ent types of scattering, as indicated in the legend. Re- 
sults obtained with both expressions MSI and MS2 of the 
form-factor are reported. The symmetry of the charac- 
teristics is related to the symmetric doping of the source 
and drain regions Q. By looking at the RS results, one 
can see that AP scattering has a negligible effect in both 
the on- and off-state regimes, while OP scattering sig- 
nificantly increases the minimum leakage current and, 
consequently, the minimum inverse subthreshold slope 
(SS), similar to the n-i-n case. As illustrated by the cur- 
rent spectrum in Fig.[7]-left, the increase of the minimum 
current is due to phonon-assisted BTBT at the source- 
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N = 13, L„ = 17 nm, V_ = 0.4 V 
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FIG. 6. Turn-on characteristics of p-i-n FETs. Top: device 
with W = 1.5 nm at Vds = 0.4 V. Bottom: device with 
W = 5 nm at Vds ~ 0.1 V. Resuhs are obtained by the 
RS and UGMS methods for different transport conditions: 
balhstic transport, AP scattering (only for the device with 
W = 1.5 nm), and both AP and OP scattering. In the MS 
case, we compare the results obtained by expressions MSI and 
MS2 of the form-factor. 
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FIG. 7. Left: current spectrum J{x, E) obtained by RS 
method for the device in Fig. [6]-top with both AP and OP 
scattering and at the bias indicated in the figure. The vertical 
line indicates the position Ls of the source-channel junction. 
The two white lines are the profile of the first conduction and 
valence subbands. Right: comparison between the current 
spectra at a; = Ls obtained by RS, UGMSl, and UGMS2 
methods. 



channel and drain-channel junctions: although the first 
conduction and valence subbands do not face each other, 
electrons can transmit from the valence to the conduc- 
tion subband by absorption of optical phonons. These 
results are in agreement with the ones in Ref. [l3l 

Again, the only difference (of about a factor of 1.5) 
between the RS and MSI results is noticed in the pres- 
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FIG. 8. Turn-on characteristics for the same device as in 
Fig. [BJ-top comparing different approximations for Sp. Top: 
only AP scattering is included, with or without the Hermitian 
part of Sp. Bottom: only OP scattering is included, with or 
without the Hermitian part of Sp. The ballistic curve is 
shown for reference in both figures and the solution method 
for each curve is indicated in the legend. 



ence of OP scattering at the point of minimum current. 
However, the accuracy with respect to RS can be almost 
completely recovered by the MS2 method, which uses 
the exact expression of the form- factor in (IT^ . This is 
more clearly shown by the comparison in Fig. [7]-right be- 
tween the current spectra at the source-channel junction 
obtained with the different methods. The lack of accu- 
racy of (PDl) in this bias condition could be related to the 
neglect of the terms F^'l^l^{i) with n ^ m, which are ac- 
tually of the same size as the terms Pj^'^^ii) included 
in (pni) . As a drawback, a slow-down of the simulation 
by about a factor of 1.5 has been measured using MS2 
compared to MSI, which can be ascribed to the increased 
computational cost of (ITS)) . 

The above considerations apply also to the turn-on 
characteristics of a wider GNR TFET in Fig. [6}-bottom. 
For this device, a lower Vds is chosen due to the lower 
band gap. In addition, we set Na = 40 (corresponding 
to ~ 5 nm), Lq = 30 nm, and a doping concentration 
in the source and drain regions of 7 ■ 10""* dopants/atom. 
4 groups of 4 modes are used in the UGMS simulations. 
It is worth noticing that the UGMS method is still accu- 
rate for the wider ribbon, indicating that the decoupling 
in separate groups is still valid, even though the subbands 
are more closely spaced than in the Na = 13 GNR. 

To evaluate the importance of the Hermitian part of 
Sp, here denoted by 5R{Sp}, we report in Fig.[S]the / vs. 
Vgs characteristics of the iVa = 13 TFET computed with 
or without 5ft{Sp}, separately for each type of scattering. 

For AP scattering (Fig. [i-top), the neglect of ^{'S^} 
leads to an underestimation of the on-state current. This 
can be understood as follows. In general, 5R{Sp} has the 
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FIG. 9. Current spectrum (left) and integral over the slab 
of the LDOS at a; = Lm (center) and a; = (right) for the 
device in Fig.[8l-top at Vas = 0.2 V and Vds = 0.4 V. Lm is 
the mid-channel position. RS method is used. The ballistic 
subband profile is shown in the inset. 
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FIG. 10. Integral over the slab of the LDOS (left), current 
spectrum (center), and average over the slab of the net elec- 
tron concentration per carbon atom (right) for the device in 
Fig. [E]-top at Vgs = 0.6 V and Vds = 0.4 V, in the presence 
of a long-range disorder potential and in the coherent limit. 
LDOS is obtained by RS, whereas the solution method for the 
other curves is indicated in the legend. Q = 3acc/4 is half the 
area of the graphene unit cell. 



B. Disordered electrostatic potential 



effect of shifting the Hamiltonian eigenvalues [2^. How- 
ever, in GNRs, the shift is of opposite sign for energies 
above and below the GNR mid-gap due to the symme- 
try of the subband structure. The result is a decrease 
of the GNR band gap, which favors BTBT, so that a 
larger current is expected when 3fi{5]p} is included in the 
simulation. Interestingly, the two expressions of 3f?{Sp} 
give almost identical values with respect to the minimum 
current. To clarify this point, we plot in Fig. [9] the cur- 
rent spectrum and the local density of states (LDOS) 
per slab at two positions along the device, for the bias 
point corresponding to the minimum current. First, it 
can be seen that the current spectrum with AP scatter- 
ing and 5ft{Sp} ^ is larger than the ballistic one: the 
reason can be ascribed to an enhanced BTBT through 
the channel region due to the band gap narrowing effect 
mentioned above, which can be appreciated from the log- 
arithmic plot of the LDOS at the mid-channel position 
in Fig. |9l-center. Secondly, by looking at Fig. [9l-right, it 
can be noticed that the peaks of the LDOS in the source 
region with K{Sp} ^ are located at the same energy 
positions as the ones of the ballistic LDOS. On the con- 
trary, the peaks of the LDOS with = are shifted 
up in energy, resulting in a tunneling current larger than 
the ballistic one and similar to the one with K{Sp} ^ 
(Fig. ini-left). The shift of the LDOS can be attributed 
to a "loss of charge" when 3?{Sp} is set to zero [30] and 
to the combined effect of the electrostatic feedback. 



For OP scattering instead, no relevant difference is ob- 
served when including 5ft{Sp}, even at high Vgs (Fig.[H]- 
bottom). 



We focus on the Na = 13 TFET. The simulations are 
performed in a non-self-consistent way, by solving the 
NEGF equations with a fixed electrostatic potential. We 
take the electrostatic potential as the sum of a disorder 
potential and the one calculated self-consistently in the 
absence of disorder and in the ballistic limit at Vgs = 
0.6 V and Yds = 0.4 V, using the RS approach. Two 
types of disorder are considered: a "long-range" one, i.e. 
slowly varying on the atomic scale, and a "short-range" 
one, i.e. rapidly varying from one atom to its neighbor 
ones. According to the first model, which is derived from 
Ref. '31', the disorder energy potential Vi at the atomic 
site i, located at position r^, is calculated as 



N 



Vi = ^S'jfjyexp I - 



(21) 



where N, I, dV are parameters, while Sj and Xj are ran- 
dom variables: Sj can take values ±1 with equal prob- 
ability; Xj is uniformly distributed over all the atomic 
position fi. The second model we study is the Ander- 
son type of disorder [s^, according to which Vi = Yi, 
where Yi is a random variable uniformly distributed in 
[-~5V/2,SV/2]. 

The coherent case, i.e. without phonon scattering, 
is considered first. In Fig. [TUI-left we show the LDOS 
corresponding to a realization of the long-range type of 
disorder, obtained with N — O.OliVc, I — 5acc, and 
SV — 0.5 eV, where Nc is the total number of carbon 
atoms inside the device and Occ is the carbon-carbon 
bond length. The resonant states induced by disorder 
are clearly visible in the LDOS. In Fig. [TUI-center and 
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FIG. 11. Same as in Fig. [TO]but in the presence of a short- 
range disorder potential. 



Fig. [TUl-right we compare the current spectrum and net 
electron charge along the device, respectively, obtained 
with RS, UGMS, and CMS. The UGMS simulations are 
performed with the same 2 groups of 4 modes used in 
the case without disorder, while the CMS ones use the 
same 8 modes but all coupled in one group. It turns out 
that the UGSM method looses accuracy due to disorder- 
induced mode mixing, while the CMS method leads to 
results very close to the RS ones (e.g. the error on the 
current value is less than 2%). 

Similar considerations can be made regarding the re- 
sults obtained with the second model of disorder with 
SV = 0.5 eV. (Fig.HH). The more regular LDOS pattern 
and the higher current spectrum in Fig. [TT] compared to 
Fig. [TU] indicate that the amount of scattering in this 
device is smaller. 

The simulations have been repeated including AP and 
OP scattering (Fig. [T2]) . For both types of disorder, it is 
seen that the MSI and MS2 approximations of the form- 
factor provide similar results and that the CMS approach 
is still accurate with respect to RS. By comparing the 
current spectra in Fig. [12] with the ones in Figs. [TOl ^ fTTl it 
can be noticed that phonon scattering adds a significant 
broadening. In the case of short-range disorder, phonon 
scattering only slightly decreases the current (from 0.83 
to 0.80 /J.A). On the other hand, in the case of long-range 
disorder, the current is increased from 0.084 to 0.10 /xA 
when phonon scattering is included in the simulation, 
indicating that the localization transport regime [29| that 
occurs in the coherent approximation is broken by the 
dephasing effect of phonon scattering. 



IV. CONCLUSIONS 

A mode space method for TB NEGF simulations of 
armchair GNR FETs including phonon scattering has 
been presented and tested with reference to both conven- 
tional and tunnel FET structures. When no disorder is 
included in the simulation, an efficient decoupling of the 
modes in different groups (UGMS) can be employed with 
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FIG. 12. Top: current spectrum at the source-channel junc- 
tion (left) and average over the slab of the net electron con- 
centration per carbon atom (right) for the device in Fig.|6]-top 
at Vas ~ 0.6 V and Vds = 0.4 V, in the presence of a long- 
range disorder potential and with both AP and OP scattering. 
The solution method for each curve is indicated in the legend. 
Bottom: same as in (top) but in the presence of a short-range 
disorder potential. 



excellent accuracy. Despite the decoupling, the method 
correctly accounts for inter-subband scattering. Simpli- 
fied expressions of the scattering-self energies have been 
compared. The one obtained by neglecting some entries 
of the form-factor is found to be accurate except for the 
bias points where transport occurs by phonon-assisted 
BTBT. On the other hand, the real part of the scatter- 
ing self-energy has only a limited effect on the device 
characteristics, especially for the case of optical phonon 
scattering, where its calculation is most demanding. In 
the presence of a disorder potential, the modes need to 
be coupled in a single group (CMS) to account for mode 
mixing, but no additional modes, compared to the ones 
used to simulate the case without disorder, need to be 
included to achieve accurate results. 

While the computational advantage of CMS over RS is 
about a constant factor with respect to the GNR width 
(equal to about 30 for the simulation parameters chosen 
in this paper), the speed-up of UGMS compared to RS 
is about 40 for a device width of 1.5 nm and increases 
proportionally to the second power of the GNR width 
(for the 5-nm-wide device considered in this paper such 
speed-up is about 360). 
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Appendix A: Calculation of the retarded phonon 
self-energy 

We consider here only the RS case, the generalization 
of the expressions to MS being straightforward. 

Replacing (U) and ([5|) into the last member of (O and 
using the identity analogous to (O valid for the G ma- 
trices, one can write 



■•AP 



■'OP 



(E) 



(Al) 



where the summation extends over all discrete energy 
points Ej e [-Emin, -Einax]- An cxprcssion of ^qp analo- 
gous to (|A3|1 containing instead of can also be 
derived. Expression (jA3l) is preferred for E higher than 
the contact Fermi energies, since the numerical error in- 
troduced by truncating the upper limit of the integral to 
-Emax is minimized due to the decaying nature of at 
high energies. On the contrary, the alternative expression 
of which depends on G> is used for low energies be- 
low the contact Fermi levels, for analogous reasons. For 
intermediate energies an average of the two formulations 
is used. 

In Sec. nil Al the results obtained with the full expres- 
sion of and S^p are compared with approximate 
solutions obtained by neglecting the respective Hermitian 
parts, i.e. 



with 



and 



i:'^p{E)^Da,JoG''{E) 



(A2) 



i:Bp{e) = Dopi o 

{Nop + l)G"iE - fiLOop) + NopG'^iE + hujop) 
G<{E - hujop) -G<{E + hujop) 



^Ue) 



Dapl O 



G>iE)-G<iE) 



and 



(A5) 



+ 00 



+iP 



G<{E' - huop) - G<{E' + hujop) 
2Tr{E - E') 



dE'} (A3) 



which does not contain G^ . It is worth noticing that the 
principal part integral in (IA3| contains only a fraction of 
the Hermitian part of Sqp [1^. It is calculated here by 
means of a piecewise constant approximation of G^(£") 
over the energy domain [-Emin , -Emax] , which is discretized 
uniformly with energy steps A e typically of the order of 
TiWop/lOO, i.e. 



OP 



(E) 



Dr. 



Io{ (Nop + l)G> (E ~ huop) 



+NopG>{E + hcuop) - [Nop + l)G<{E -f fiujo 



^NopG<iE-huJop) 



(A6) 



+ 00 



G<{E't^lUo 



E-E' 



-dE' 



Y.G<{E,)\n 



E-E^+AE/2Thio. 



op 



E -Ej- Ab/2=f ?iWo 



Expression (|A6[) is much more efficient than 
(A4) the absence of the integral term. 



31) due to 



[1] M. Y. Han, B. Ozyilmaz, Y. Zhang, and P. Kim, Phys. 
Rev. Lett. 98, 206805 (2007). 

[2] Z. Chen, Y.-M. Lin, M. J. Rooks, and P. Avouris, Phys- 
ica E 40, 228 (2007). 

[3] X. Wang, Y. Ouyang, X. Li, H. Wang, J. Guo, and 
H. Dai, Phys. Rev. Lett. 100, 206803 (2008). 

[4] G. Fiori and G. lannaccone, IEEE Electron Device Let- 
ters 28, 760 (2007). 

[5] Y. Ouyang, Y. Yoon, and J. Guo, IEEE Trans. Electron 
Devices 54, 2223 (2007). 



[6] G. Liang, N. Neophytou, M. S. Lundstrom, and D. E. 

Nikonov, J. Appl. Phys. 102, 054307 (2007). 
[7] Y. Yoon, G. Fiori, S. Hong, G. lannaccone, and J. Guo, 

IEEE Trans. Electron Devices 55, 2314 (2008). 
[8] P. Zhao, J. Chauhan, and J. Guo, Nano Lett. 9, 684 

(2009). 

[9] R. Grassi, A. Gnudi, E. Gnani, S. Reggiani, and G. Bac- 
carani, J. Comput. Electronics 8, 441 (2009). 
[10] S. Datta, Quantum Transport: Atom to Transistor 
(Cambridge University Press, Cambridge, UK, 2005). 



[11] Y. Ouyang, X. Wang, H. Dai, and J. Guo, Appl. Phys. 

Lett. 92, 243124 (2008). 
[12] Y. Yoon, D. E. Nikonov, and S. Salahuddin, Appl. Phys. 

Lett. 98, 203503 (2011). 
[13] Y. Yoon and S. Salahuddin, Appl. Phys. Lett. 101, 

263501 (2012). 

[14] N. D. Akhavan, G. Jolley, G. A. Umana-Membreno, 
J. Antoszewski, and L. Faraone, J. Appl. Phys. 112, 
094505 (2012). 

[15] J. Wang, E. Polizzi, and M. Lundstrom, J. Appl. Phys. 

96, 2192 (2004). 
[16] M. Luisier, A. Schenk, and W. Fichtner, J. Appl. Phys. 

100, 043713 (2006). 
[17] S. Jin, Y. J. Park, and H. S. Min, J. Appl. Phys. 99, 

123719 (2006). 

[18] S. Poll, Modelling and simulations of postCMOS devices, 
Ph.D. thesis. University of Bologna (2009). 

[19] A. Esposito, M. Frey, and A. Schenk, J. Comput. Elec- 
tron. 8, 336 (2009). 

[20] D. Nikonov, H. Pal, and G. Bouri- 

anoff, "Scattering in NEGF: Made simple," 
|http ■ //na nohub . org/resources/7772 (2009) . 

[21] A. Afzalian, J. Appl. Phys . Ti0TiO70945 1 7 (2011). 



10 



[22] M. Shin, J. Appl. Phys. 106, 054505 (2009). 

[23] R. Grassi, A. Gnudi, E. Gnani, S. Reggiani, and G. Bac- 

carani, IEEE Transactions on Nanotechnology 10, 371 

(2011). 

[24] P. Zhao and J. Guo, J. Appl. Phys. 105, 034503 (2009). 
[25] Y.-W. Son, M. L. Cohen, and S. G. Louie, Phys. Rev. 

Lett. 97, 216803 (2006). 
[26] R. Lake, G. Klimeck, R. C. Bowen, and D. Jovanovic, J. 

Appl. Phys. 81, 7845 (1997). 
[27] A. Svizhenko, M. P. Anantram, T. R. Govindan, 

B. Biegel, and R. Venugopal, J. Appl. Phys. 91, 2343 

(2002). 

[28] S. O. Koswatta, M. S. Lundstrom, M. P. Anantram, and 
D. E. Nikonov, Appl. Phys. Lett. 87, 253107 (2005). 

[29] S. Datta, Electronic Transport in Mesoscopic Systems 
(Cambridge University Press, Cambridge, UK, 1997). 

[30] A. Svizhenko and M. P. Anantram, Phys. Rev. B 72, 
085430 (2005). 

[31] M. Poljak, E. B. Song, M. Wang, T. Suligoj, and K. L. 

Wang, IEEE Trans. Electron Devices 59, 3231 (2012). 
[32] A. Lherbier, B. Biel, Y.-M. Niquet, and S. Roche, Phys. 

Rev. Lett. 100, 036803 (2008). 



